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We investigate the rotational dynamics of disk accretion around a strongly 
magnetized neutron star with an aligned dipole field. The magnetospheric field 
is assumed to thread the disk plasma both inside and outside the corotation ra- 
dius. As a result of disk-star interaction, the magnetic torque on the disk affects 
Oh" the structure of accretion flow to yield the observed spin-up or spin-down rates 

for a source of given fastness, magnetic field strength, and mass accretion rate. 
Within the model we obtain a prescription for the dynamical viscosity of such 
magnetically modified solutions for a Keplerian disk. We then use this prescrip- 
tion to find a model solution for the rotation rate profile throughout the entire 
^ ■ disk including the non-Keplerian inner disk. We find that the non-Keplerian 

angular velocity transition region is not necessarily narrow for a source of given 
spin state. The boundary layer approximation, as in the standard magnetically 
threaded disk model, holds only in the case of dynamical viscosity decreasing all 
the way to the innermost edge of the disk. These results are applied to several 
observed disk-fed X-ray pulsars that have exhibited quasi-periodic oscillations 
(QPOs). The QPO frequencies provide a constraint on the fastness parameter 
and enable one to determine uniquely the width of the angular velocity transition 
zone for each source within model assumptions. We discuss the implications of 
these results on the value of the critical fastness parameter for a magnetized star 
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in spin equilibrium. Applications of our model are also made with relevant pa- 
rameters from recent numerical simulations of quasi-stationary disk-magnetized 
star interactions. 

Subject headings: accretion, accretion disks — magnetic fields — stars: magnetic 
fields — stars: neutron — X-rays: stars 

1. INTRODUCTION 

The evidence for the interaction between a rotating magnetized star and a surrounding 
accretion disk emerges in a wide variety of astrophysical systems such as X-ray binary pulsars 
(Joss & Rappaport 1984), magnetic white dwarfs in cataclysmic variables (Warner 1995), 
and T Tauri stars (Basri & Bertout 1989). The detailed modelling of this interaction would 
contribute to our understanding of many problems like the time rate of change in the pulse 
periods of the disk fed neutron stars and the quasi-periodic oscillations (QPOs) observed in 
some sources (Alpar & Shaham 1985; Aly & Kuijpers 1990). 

In the early investigations (Scharlemann 1978; Ghosh & Lamb 1978; Aly 1980), the 
degree of diamagnetism presumed for the disk plasma has led to very different magnetic field 
configurations. If the disk is fully diamagnetic, then it excludes the stellar field completely 
and the accretion flow from the disk midplane to the surface of the star is assumed to start 
within a very narrow region near the inner disk radius (Ichimaru 1978). If the accreting 
plasma has a non-zero resistivity, then the stellar field may penetrate the disk and thread it 
in a broad zone via the Kelvin- Helmholtz instability, magnetic field reconnection with small 
scale fields in the disk and turbulent diffusion (Ghosh & Lamb 1979, hereafter GL79; Wang 
1987, hereafter W87). One of the most important predictions of the magnetically threaded 
disk (MTD) model is that the star can spin down while matter accretes on it. 

In the MTD model of GL79, the star-disk interaction region consists of two distinct 
parts. In a broad outer transition zone, where the angular velocity is Keplerian, the effective 
viscous stress is dominant compared to magnetic stress associated with the twisted field 
lines. The disk matter is brought into corotation with the neutron star only in a narrow 
inner transition zone or boundary layer within which magnetic stress dominates over viscous 
stress. The disk resistivity is very high so that the magnetic diffusivity exceeds the turbulent 
viscosity by several orders of magnitude. The coronal plasma outside the disk is force-free 
and the disk is threaded over a large range of radii by magnetic field lines that close through 
the neutron star. 

The study of field line twisting in a force-free magnetosphere has revealed that a closed 
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field line tends to innate and evolve into an open one breaking the disk-star link (Lynden- 
Bell & Boily 1994; Lovelace, Romanova, & Bisnovatyi-Kogan 1995, hereafter LRBK95). 
The strong outflows associated with the opening of magnetic field lines were reported in 
several two dimensional magnetohydrodynamic (MHD) simulations (Hayashi, Shibata, & 
Matsumoto 1996; Miller & Stone 1997; Goodson, Winglee, & Bohm 1997; Goodson, Bohm, 
& Winglee 1999). Although these numerical studies indicated that the accretion flow around 
a strongly magnetized object can be nonstationary and episodic in character, the recent MHD 
simulations by Romanova et al. (2002, hereafter RUKL02) showed that quiescent accretion 
regimes are also possible provided a matter dominated differentially rotating corona instead 
of a magnetically dominated force-free one is used as an equilibrium initial condition to 
examine quasi-stationary situations. The main reason for obtaining different results in the 
above mentioned simulations is that in the case of a corona corotating with the central star (or 
non- rotating corona), the initial discontinuity in the angular velocity between the Keplerian 
disk and the magnetosphere leads to huge field twisting on the corona-disk boundary. This 
gives rise to the generation of very large toroidal field and to strong magnetic braking of the 
disk. A differentially rotating matter dominated corona is an appropriate initial condition to 
avoid this discontinuity (Romanova et al. 1998). It reduces the initial magnetic braking of 
the disk and allows one to investigate the quiescent phase of disk-star interaction in long-term 
numerical simulations (RUKL02). 

The physical mechanism responsible for the disruption of the inner region of an axisym- 
metric disk around a magnetic rotator was investigated by a number of authors (Brandenburg 
& Campbell 1998; Campbell & Heptinstall 1998, hereafter CH98). Integrating the relevant 
MHD equations radially from the weakly magnetic outer disk regions inwards, they found a 
vertically disrupted and viscously unstable disk solution inside a critical radius. They also 
argued that the viscous instability associated with the elevated temperature in the inner 
disk regions is due to magnetically enhanced viscous dissipation, which causes the vertical 
equilibrium to break down when the radiation pressure becomes significant. We point out in 
the present paper (see §§ 3 and 4), that the viscous instability proposed by these authors as 
the mechanism responsible for the truncation of the inner disk is a consequence of choosing 
Keplerian rotation for the disk plasma throughout the accretion flow, artificial in the sense 
that the angular motion could easily be modified by the huge viscous stress estimated for 
the inner parts of the disk. 

For a disk around a non- magnetized star, Glatzel (1992) solved this problem without 
making the assumption that the rotation rate is Keplerian except in a very narrow boundary 
layer. In this paper we will present the analogous solution for the magnetic case, allowing for 
non- Keplerian rotation in the inner disk. Numerical simulations by RUKL02 also have shown 
non-Keplerian rotation in the magnetically braked inner regions of the disk. We consider 
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an axisymmetric thin disk threaded by the magnetic dipole field of an aligned rotator and 
explore the effect of the rotating magnetosphere on the angular velocity profile of the disk 
plasma. In § 2 we introduce the relevant MHD equations to derive a dimensionless function 
for the vertically integrated dynamical viscosity. In § 3 we discuss the behavior of this 
function in connection with the net torque acting on the central star. In § 4 we employ 
the dynamical viscosity and solve the equation of angular momentum conservation for the 
rotation profile of a disk around a neutron star in spin equilibrium without making boundary 
layer approximation. We apply our model to several X-ray binary pulsars that have exhibited 
QPOs. Using the observed periods, X-ray luminosities, QPO frequencies, and spin-up or 
spin-down rates of these sources, we obtain constraints on the model parameters within 
the beat frequency model and estimate the appropriate disk rotation curve for each source. 
We find that the inner disk rotation rate adapts to Keplerian rotation in the outer disk 
through a transition zone whose width and rotation rate profile are such as to accomodate 
the mass accretion rate, rotation period, and magnetic moment of the central object and 
to yield the observed spin-up or spin-down rates. In § 5 we concentrate on the qualitative 
similarities between our model and the results of numerical simulations by RUKL02 on the 
disk structure. Using their numerical data, we roughly estimate the radius where the disk 
structure is significantly changed and magnetic braking is efficient. We summarize and 
discuss our results in § 6. 



2. BASIC EQUATIONS 

The steady-state equation of conservation of angular momentum for an axisymmetric 
magnetized disk around a neutron star having a dipole moment aligned with its rotation 
axis can be written in cylindrical coordinates (r, 0, z) as 
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where v is the coefficient of viscosity. The conservation of mass for a steady-state, axisym- 
metric disk is given by the continuity equation, 
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which can be vertically integrated to yield the constant mass influx condition, 

-2tt I rpv r dz = M , (3) 

J-H 

provided that pv z is negligible at z = ±H, a valid assumption except for the case of strong 
winds or magnetically driven outflows. Here, M is the constant mass inflow rate and H is 
the half-thickness of the disk. 

In a geometrically thin accretion disk (i.e., H Cr), the angular momentum balance (1) 
can be approximated as 

— (rpv r r — [rpv z r Q) = — \vpr — j 

+ 4vr dz ' K V 

where = fir is used and the r<p component of the magnetic stress is neglected in comparison 
with the last term in equation (1), for the ratio of their magnitudes in the thin disk limit 
gives 

\rB r d r (rBs)\ 



B r 



B 7 
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provided \B r /B z \ r/H, a condition that can only be violated for extremely large disk 
conductivities. 

The vertically averaged form of the angular momentum equation (4), together with the 
use of equation (3), is finally obtained as 

Tr (^^f + A:fr2 ") = ~ rlB - B i • (6) 

where B^(r) = B^r, z — H) — —B<p(r, z = —H) due to the field antisymmetry and E is the 
surface density function defined by 

S(r) = f p(r,z)dz = 2pH , (7) 

J-H 

where p denotes simply the vertically averaged density. The last term in equation (1) was 
also identified as the dominant magnetic term in the early derivation of equation (6) by 
Campbell (1992) for steady disks and by Lovelace, Romanova, & Newman (1994) in a form 
appropriate for time-dependent disks with outbursts. 

The vertical component of the magnetic field threading the disk can be expressed, in 
general, as 

Bz (r) = - S (r)/i,r" 3 , (8) 
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where //* = (1/2) B*Rl is the dipole moment of the neutron star of radius i?* in terms of 
its polar surface field strength B* and s(r) represents a screening coefficient which accounts 
for the effect of induced currents on the pure stellar dipole field in a partially diamagnetic 
disk. In a self-consistent treatment, the radial dependence of s must be properly found from 
a detailed analysis of the magnetospheric current system in the corona. The modelling of 
such a disk-corona-star system is beyond the scope of the present study. In this work, we 
will rather assume that the vertical component of the magnetic field threading the disk can 
be represented by 



B z (r) = -s cff 



At* 



(9) 



where s e ff is an effective screening coefficient taken to be a constant over a wide range of 
radii (Wang 1995). 

The azimuthal component of the magnetic field at the surface of a disk in vertical 
hydrostatic equilibrium (i.e., v z = 0) can be estimated from the toroidal component of 
the induction equation, which can be written under the assumptions of steady-state and 
axisymmetry as 

-(fir^) + - (nrB r ) - - (vrB,) 
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where rj is the magnetic diffusivity of the disk. Assuming that all quantities change on a 
length scale r in the radial direction and H in the vertical direction, the significance of the 
second term compared with the first term in equation (10) can be estimated as 
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provided we have \B. 
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provided r\ > \v r \ H is satisfied. The consistency of equation (11) with equation (13) can be 
justified using the poloidal component of the induction equation which can be written as 
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where the disk is assumed to be in vertical hydrostatic equilibrium, i.e., v z = 0. For very 
small radial inflow velocities, i.e., v r <C i]/r, equation (14) yields \B r \ ~ (H/r) \B Z \ and the 
distortion of the external field due to accretion becomes negligible. In this case the estimates 
in equations (11) and (12) are both 0((H/r) 2 ) <C 1, while the estimate in (13) is smaller 
than 0((H/r) 2 ) by the factor v r rjr\. If, on the other hand, \B r \ < |-B 2 |, as is likely when the 
field is distorted, then \d r B z \ <C \d z B r \ and we obtain 

Bt = ~—B z 15 
V 

from the vertical integration of equation (14), where B+(r) = B r (r, z — H) — —B r (r,z = 
—H) following the field antisymmetry (see Lovelace, Romanova, & Newman 1994; also 
LRBK95). When equation (15) holds, approximations (11) and (13) are both < 0(H/r) <C I. 

Taking the conditions (11), (12), and (13) to be satisfied, the toroidal component of the 
induction equation (10) can be safely reduced to 

„ on of DBA , x 

^a:+fl:iir =°- ( 16 ) 



dz dz \ dz 

The vertical gradient of angular velocity depends on the scale height over which the plasma 
is brought into corotation with the neutron star. In the case of a force-free magnetosphere 
outside the disk, this scale height can be estimated as H, the half-thickness of the disk, 
since the magnetic energy density becomes much greater than the kinetic energy density 
of the plasma for \z\ > H. The resulting twisting of field loops due to this huge vertical 
shear gradient may lead to an open field line configuration unless \d£L/dz\ < — Q\/H. 
The recent time- dependent simulations by RUKL02 have shown that the field line opening 
is strongly suppressed for a relatively dense corona outside the disk. Assuming that the 
vertical shear along a flux tube linking the star and the disk is reduced by the inertial effect 
of a corona, the integration of equation (16) over z can be approximated as 

B+ = er (0, - n) Hi]- X B Z , (17) 

where the transition from the angular velocity of the disk plasma, Q, to the rotational rate 
of the neutron star, Q*, takes place within an effective scale height of t~ x H with e < 1 (see 
Campbell 1992; also LRBK95 for other derivations of eq. [17]). Here, e is a shear reduction 
factor which may arise from a differentially rotating corona to result in a small field twist 
at the surface of the disk. Even if the corona is nearly force-free, it cannot be current-free. 
As the magnetic torque is transmitted to the neutron star by the poloidal currents flowing 
across the disk surface through the magnetosphere, the rigid-body rotation of the corona 
may not be realized (see W87; also CH98 for a similar shear reduction factor). 
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We consider a specific disk region (e.g., a narrow boundary layer of width 5r < H <C r) 
where the contribution from the radial gradient terms in equation (10) is not negligible. In 
this case, the angular velocity transition from Q to occurs within 5r. The approximate 
balance between the second term (field generating) and the third (advection) and fourth 
(dissipation) terms in equation (10) yields a crude estimation for the toroidal field component, 

B+ « r (a - Q) (5r) rf x B+ ~ r (0, - Q) Erf x B z , (18) 

provided the width of the boundary layer is assumed to be the electromagnetic screening 
length for the poloidal magnetic field, i.e. Sr ~ rj/v r (see GL79). The vertical gradient terms 
are not taken into account for simplicity in the radial integration of equation (10). The last 
step in (18) follows from / B z ~ H/5r according to equation (14). This analysis indicates 
that the form of the expression for the azimuthal field component in equation (17) may 
generally be applicable for the entire disk regardless of the radial extension of the angular 
velocity transition zone. 

In the thin disk limit, the vertically integrated forms of the r and z components of the 
momentum equation can be written as 

(nl-tf)rH-^^O, (19) 

where c s = (P/ p) 1 ^ 2 is the speed of sound in terms of the mid-plane temperature of the disk, 
and fix is the Keplerian angular velocity given by 

n K (r) = (^J , (21) 

where M* is the mass of the neutron star. The magnetic tension force in equation (19) 
is the main agent responsible for the deviation of Q(r) from f2x(r) if the accretion disk is 
thin and the poloidal magnetic field varies on a length scale r. The vertical equilibrium 
of a magnetically threaded thin disk is given by equation (20), where the gas pressure is 
balanced by the magnetic pressure of the horizontal field components in addition to the 
vertical component of gravity. Using equations (19) and (20), it is possible to make an order 
of magnitude estimate for the sound speed. In the absence of a strong global magnetic field, 
the rotational profile is nearly Keplerian and c s ~ £l K H. The rotation of the disk matter at 
a sub-Keplerian rate can only be realized provided we have 

J^i£ < m Hr (22) 
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in equation (19) with \B Z \ ~ B. If the magnetic diffusivity is turbulent in nature, then the 
radial component of the magnetic field can be significant, i.e., \B+\ < B, and the vertical 
equilibrium (20) is satisfied by the balance between the gas and magnetic pressures alone, 

C> & P ~ ^ ' (23) 
since the gravity term is relatively small in a thin disk, i.e., Q^H 2 <C Q^Hr. 

After equation (23), we propose the following prescription for the sound speed: 

c 2 s = £(r) tf K Hr , (24) 

where £(r) < 1 is a dimensionless factor of unknown radial dependence. Its explicit form can 
only be deduced from a detailed analysis of the field configuration and density distribution 
in the disk. 

In this work, we assume that the magnetic diffusivity of the disk is of the same physical 
origin (i.e., turbulence) as the viscosity, and we adopt the original a formalism (Shakura & 
Sunyaev 1973, hereafter SS) to write 

V = ad-^, (25) 

where a rf < 1 is a dimensionless numerical factor. 

Combining equations (24) and (25), the azimuthal field component at the disk surface 
(see eq. [17]) is simplified to 

b+ = l4> (a, - n) n^B z , (26) 

where 7^ = e/^ctd will be treated as a constant for simplicity (see, e.g., Livio & Pringle 1992). 
As we will see in § 3, the vertically integrated dynamical viscosity, z/£, for a magnetically 
threaded disk changes with distance. The radial variation of uT, can be achieved in general 
provided both the coefficient of viscosity v and the surface density £ are functions of radius. 
As the viscosity v and diffusivity rj have the same, turbulent nature by assumption, the 
magnetic Prandtl number, v/i], should be of order one at different radii of the disk. This also 
requires that the parameter £ should have a radial dependence in general for self-consistency. 
In this case, 7^ would not be approximated as a constant unless the shear reduction factor 
e changes with distance in more or less the same way as £ does. Our physical motivation 
for taking 7^ as a constant is the following. The compressive effect of the magnetic field 
on the vertical equilibrium of the disk becomes negligible in the weakly magnetized outer 
disk regions. The gas pressure is mainly balanced by the gravitational force and c s ~ Q K H 



-10- 



(see eq. [20]). According to equation (24), this is equivalent to choose £(r) if/r in the 
outer disk. The resulting toroidal field is 0(r / H)B Z ^> B z unless e(r) ~ H/r. The shear 
reduction factor e is therefore expected to decrease from 0(1) at small radii to 0(H/r) at 
large radii provided the disk is threaded by the stellar field lines that are closed and stable 
at least to some extent beyond the corotation radius. As there appears to be a correlation 
between e(r) and £(r), we approximate their ratio as a constant throughout the entire disk. 

In the absence of a complete theory of the disk-star magnetic interaction, we will not 
attempt to obtain a full disk solution at this stage; we will rather concentrate on the con- 
servation of angular momentum to treat the dynamics of disk accretion under the action of 
an external magnetic field of stellar origin. Using equation (26), it follows from equation (6) 
that 

^ (27rz/Sr 3 ^ + Mr 2 ^ = - l4> r 2 (Q* - Q) tl^B 2 (27) 

which can be solved for the rotation profile, £l(r), for a given vertically integrated dynamical 
viscosity, z/E, if B z (r) is known. 



3. DYNAMICAL VISCOSITY 

In the following, we will deal with the angular momentum balance (see § 2) to derive 
an expression for the vertically integrated dynamical viscosity. Before attempting to solve 
equation (27), we scale the variable quantities r, z/X, and Q by their typical values such 
that x = r/r in is a dimensionless coordinate in units of the inner disk radius r in where we 
assume f2 = f(x) = 3ttz/E/M is a dimensionless function for the dynamical viscosity, 
and uj{x) = f2/f2K(r in ) is a dimensionless angular velocity for the disk plasma. Finally, we 
use equation (9) for the screened dipole field to rewrite equation (27) in a non-dimensional 
form: 

a ro a, , 

-P K - lu) x^ 12 , (28) 



d_ 

dx 



2 „, x 3 dw 2 
-f(x)x- + xu 



where cu* = Jl*/ilK(?"in) = (?"coAm) 3 ^ 2 with the corotation radius r co = (GM,/^) 1 / 3 and 
(3 = 70 s 2 fi(rx/r- n) ) 7 l 2 with the Alfven radius given by 



r A = M- 2 ' 7 ^J 7 {GM, 



-1/7 



-1/7 



3.4 x 10" cm M^^til ( ) • (29) 



Here, Mn is the mass accretion rate expressed in units of 10 17 g s 1 and /^3o is the magnetic 
dipole moment in units of 10 30 G cm 3 . If the neutron star magnetic field is sufficiently weak 
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for a given mass accretion rate, e.g., £?* < 10 7 G in equation (29), then we have r A < R* 
and the disk may extend down to the surface of the central object (i.e., r in = R*). For such 
weakly magnetized systems also known as the standard a-disks (SS), we get (3 — 0. The same 
is true also if the disk plasma is perfectly diamagnetic. In the conventional picture of a-disks, 
the rotation curve is Keplerian throughout the accretion flow except in a narrow boundary 
layer situated at the inner edge of the disk where the angular velocity of the plasma changes 
from f2 K (-R*) to Q*. In a standard thin disk around a non-magnetic star, we can integrate 
equation (28) with (3 = to find f(x) = 1 — (x^ 1 / 2 using the Keplerian profile, oo(x) = x^l 2 . 
The integration constant (, here, denotes the inflow rate of angular momentum J, across a 
cylindrical boundary at radius r, in terms of the angular momentum flux carried by matter 
onto the neutron star through the inner edge of the disk, i.e., J = (MRIQk(R*). According 
to the standard approach (Lynden-Bell & Pringle 1974), there is no viscous torque acting 
on the inner disk. This corresponds to ( ~ 1 and the boundary layer is extremely thin 
in extension. However, the angular velocity transition region is not necessarily narrow as 
Glatzel showed (Glatzel 1992, hereafter G92) if ( is treated as a free parameter in f(x), which 
in turn can be used to solve for the structure of the accretion flow (Fujimoto 1995). In the 
present investigation, we will derive a form of vertically integrated dynamical viscosity in a 
way similar to the work of G92 but generalized to treat the disk interacting with a magnetic 
star. In our case, /3^0 and the angular momentum exchange between the magnetosphere 
and the disk is not negligible except for sufficiently large distances from the magnetically 
braked inner disk regions. The Keplerian rotation can only be regarded as a particular 
solution that can be employed to obtain a prescription for the dynamical viscosity, which 
may hold throughout the accretion flow provided the angular motion of the inner disk plasma 
can be matched to the rotation rate of the stellar magnetosphere. The radial integration of 
equation (28) for a Keplerian angular velocity profile yields 

f(x;(3,j,u;*) = 1 - ^ Wi x- 2 + l -Px- 7 ' 2 - jx^ 2 , (30) 

where j = J/ Mr j^K^m) is the net angular momentum flux into the neutron star unless 
there are some other efficient mechanisms of angular momentum loss for the disk material 
such as winds associated with the opening of field lines. Although j, here, appears to be an 
arbitrary integration constant, it is the dimensionless torque applied by the disk on the star, 
i.e., the star spins up (or spins down) for j > (or j < 0). 

A vertically integrated dynamical viscosity form similar to the one we consider in 
equation (30) was derived and used previously by Brandenburg & Campbell (1998, here- 
after BC98) to solve for the disk structure around strongly magnetic accretors. These au- 
thors, however, assumed j = to ensure that the disk structure matches that of weakly 
magnetized accretion at a sufficiently large radius r out . Adopting z/£ = M/3n, that is, 
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f(x — r out /r in ) = 1 as an arbitrary outer boundary condition, they integrated the relevant 
MHD equations in the radial direction from the outer disk regions inwards. Note that / = 1 
at large distances whether j = or not. However, the structure at finite x depends on the 
relative importance of the terms, and the value of j will be an important parameter in de- 
termining the structure. BC98 did not attempt to fit their solution to a magnetically braked 
inner accretion flow; they assumed that the Keplerian rotation holds throughout the disk. 
Their model implies that the torque on the star is always zero, i.e., j = 0, and the angular 
velocity of the plasma still remains Keplerian in the inner parts of the disk even though the 
viscous dissipation there due to magnetically enhanced viscous stress (see the third term 
on the r.h.s. of eq. [30]) leads to the divergence of the disk height as a result of elevated 
temperature. Both assumptions are very restrictive and non-physical in the sense that the 
choice of j becomes quite significant at small radii although it seems to have no effect on 
the disk structure at large radii (see § 4) and the huge viscous dissipation estimated for 
the inner disk regions results from using overestimated angular frequencies (e.g., Keplerian 
frequencies) together with a magnetically enhanced dynamical viscosity. The last argument 
can be made clear if we consider the conservation of energy for a radiatively efficient disk, 

2aT * w f H Qvdz = 1^ ( X W 2 tt x >P>i> w J ' (31) 

where T s is the effective surface temperature of the disk and Q v = vp(rd r VL) 2 is the viscous 
dissipation rate per unit volume. The approximation sign in equation (31) indicates that 
there may exist additional sources of heat (e.g., Ohmic dissipation) other than viscosity in 
the disk. If the Keplerian assumption, i.e., uj{x) = x~ 3//2 , were valid throughout the entire 
disk as in the work of BC98, then we would certainly find huge temperatures as x — > 1 for 
sufficiently high (3 and low values (see, e.g., CH98 for the corresponding values) in the 
case of j = (see also, eq. [30]). 

What we expect to find in steady accretion is that the inner boundary condition, that 
is, the fastness of the stellar magnetosphere cj*, adjusts the rotational disk dynamics for a 
given torque j. The quantity (5 in equation (30), which reflects the magnetic configuration 
in the inner disk, is not an independent parameter. Rather, it is linked intimately to j and 
its value can be calculated (see § 4) when the angular velocity profile is determined from 
equation (28). 

The definition of (3 with r in = u^ 3 r co , yields for a typical X-ray binary pulsar: 

P c 12 7 , sl s uj^MJP^Iso ( 7^£) " 5/3 > ( 32 ) 
where P* is the spin period of the neutron star. Using the torque expression, AT* = Pfi* = 
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jM(GM*r in ) 1/2 with a moment of inertia h = (2/5) M*Rl, we find 




(33) 



where P* 6 is the neutron star radius expressed in units of 10 6 cm and P* is the time rate of 



change of the spin period. 

We now examine some of the basic features of the dimensionless function we introduced 
for the dynamical viscosity (see eq. [30]) before we present various types of rotation profiles 
in the next section, for the behavior of uj{x) strongly depends on f(x) and its first derivative 
f'(x) (see eq. [28]). As an illustrative example of how the dimensionless viscosity is affected 
by the last term in equation (30), we display in Figure 1 plots of f(x) for two different values 
of j. The dynamical viscosity function shown by the solid curve in Figure 1 corresponds 
to (3 = 16, cj* = 0.6, and j = —0.8. According to equations (32) and (33 ), this can be 
realized for a typical X-ray binary pulsar of spin period 1 s if 7^ ~ 0.4 with s c g ~ 1 and 
P* ~ 1.7 x 10- 12 ss _1 . The dashed curve (see Fig. 1) with (3 = 13, = 0.3, and j = 0.9, 
represents another physically plausible dynamical viscosity for the same pulsar provided 
70 ~ 0.1 for s cS ~ 0.8 and P* ~ —1.5 x 10 _12 ss _1 . Note that the curve associated with 
a negative value of j (spin-down) is characterized by two local extrema, corresponding to 
a maximum and a minimum, located respectively at X2 and x±, with x\ < Xi- Also, note 
that there is no local maximum for j > 0. The curves in Figure 1 are plotted for illustrative 
values of (3 and cj*. The conclusions drawn from Figure 1 however can be generalized also 
for other values of (3 and a;*, if we consider the sign of the second derivative of f(x) at local 
extrema. The real extremum points, x ext = £1,2, can be found from f'(x cxt ) = as 



provided we have 16/3u; 2 + 21 j > 0. The first root, x±, specified by the plus sign within the 
parenthesis in equation (34) can be identified to be the local minimum, whereas the second 
root, X2, with minus sign represents the local maximum. Equation (34) then, implies that 
x\ < X2- The local minima of both curves in Figure 1 are inside the corotation radii, i.e., 
x 1 ~ 1.3 < x co ~ 1.4 < x 2 ~ 9.8 for uo* = 0.6 and X\ ~ 1.8 < x co ~ 2.2 for uj* = 0.3, where 
x co = r co /r in . As we will see in § 4, x\ < r co /r in is always satisfied if the rotation of the disk 
matter is nearly Keplerian for r > r co . The plasma is progressively brought into corotation 
with the neutron star in an inner transition zone of radial extension 5r- in . The angular 
velocity of the inner disk matter starts deviating from its Keplerian value at r < xir in , for 



xi,2 = ^ [(-36/3w. ± 9^16/^2 + 21/?/) j 




(34) 
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the magnetic stress dominates over the shear stress when the dynamical viscosity becomes 
sufficiently small, that is, when f(xi) <C 1. The width of the inner transition zone can be 
then estimated as S < x-y — 1. 



4. APPLICATION TO OBSERVED BINARY X-RAY PULSARS 

We solved the second order differential equation (28) for u(x), substituting equation 
(30) for the dynamical viscosity. As we do not know a priori where the true inner disk 
radius is located, we tried different values of uj* for a given torque j. As an initial guess for 
each set of (a>*, j), we calculated (3 such that f(xi) = 0. We observed that the necessary 
boundary conditions, i.e., u(x = 1) = u;* and uo(x > x co ) ~ x~ 3//2 , are satisfied to the desired 
accuracy in our numerical iteration only for certain values of f3 for which the dynamical 
viscosity nearly vanishes at x — x\. Our results for a rotator in spin equilibrium, i.e., j = 0, 
are shown in Figures 2a and 26. As in Figure 1, the corresponding dynamical viscosities 
(dashed curves) and the angular velocity profiles (solid curves) are plotted as functions of 
x = r/r- m . The dotted curves represent Keplerian rotation. Note that the angular velocity is 
braked by the magnetic stress only inside x co , and the rotation of the disk plasma is almost 
Keplerian outside x co . We calculate the net torque applied by the disk on the star using a 
closed surface through the corotation radius. The surface encloses the central object and 
excludes part of the disk for x > x co . This choice avoids the uncertainties involved in making 
torque calculations in the inner disk where Q ^ Qk- We write the torque as 

N* = + N£ + ^mag , (35) 

where the angular momentum flux carried by the material stress through a cylindrical surface 
area 4nr co H(r co ) is 

iV c r = Mr 2 co n K (r co ) . (36) 

The contribution of the viscous stress to the flux of angular momentum through the same 
surface is given by 

AC = |^o)r c 2 (^) (27rr co ) 

^ ' ^co 

= -f(x co )Mr 2 co n K (r co ) . (37) 
The magnetic coupling outside the corotation radius yields a net spindown torque, 

poo 1 

AW = - / r 2 B z B+ dr = --fe 7 / 2 Mr c 2 o K (rc„) (38) 
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which, together with equations (36) and (37), gives iV* = jM(GM*r hl ) l l 2 according to 



equation (30). This justifies our interpretation of j as the dimensionless torque acting on 
the neutron star in units of the angular momentum flux carried by matter through the inner 
edge of the disk. Also, note from Figures 2a and 2b that x\ decreases (i.e., X\ — > 1) while 
u;* increases as expected. There is, however, an upper limit for if the star is in rotational 
equilibrium since 



independent of /3, following from equation (34). Thus, for accreting systems in rotational 
equilibrium, j = 0, the critical fastness parameter u;* jC for the transition from spin-down to 
spin-up depends on the angular velocity profile (width of the transition zone) in equilibrium. 
Equation (39) implies that the maximum possible critical fastness parameter is (c^*, c )max = 
0.875, a value quoted previously by Wang (1995, hereafter W95) as the critical fastness 
parameter for a star accreting from a turbulent disk. According to the accretion torque 
model of W95, the Keplerian rotation holds throughout the whole disk. Equivalently, in our 
model, this maximum critical fastness parameter obtains if the angular velocity transition 
region is extremely narrow or there is no transition region at all. This corresponds to the 
x± — 1 limit among systems in spin equilibrium as seen from the comparison of Figures 2a 
and 2b. The usual fastness parameter, u s , can be expressed in terms of cj* as 



Here, tq is the radius where Q reaches its maximum and xq = ro/r- m . Using equation (40), 
we obtain u; SjC = 0.75 for o;*^ = 0.5 from Figure 2a and cj sc = 0.83 for ^ c = 0.8 from 
Figure 2b. Our analysis indicates broad transition zones for j = rather than narrow ones 
if cj SjC < 0.8. The disk is Keplerian at all radii only for o;* iC = u SiC = 0.875 and xo — x\ — 1. 

We have extended our calculations to several binary X-ray pulsars observed as QPO 
sources. One of the simplest and most frequently used models for QPOs from X-ray pulsars 
is the beat frequency model (Alpar & Shaham 1985). According to the beat frequency 
model (BFM), the inhomogeneities at the inner edge of the Keplerian disk modulate the 
X-ray intensity and the QPO is observed as the beat frequency between the local Keplerian 
frequency and the neutron star spin frequency, 




(39) 




(40) 



^QPO = M r o) - ^* • 



(41) 



The corresponding fastness parameter can be written as 

, bfm 1 



(42) 



1 + Z/qpo ' 
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where P* is the spin period of the observed X-ray pulsar. The QPOs detected from X-ray 
pulsars can be used, within the beat frequency model, to obtain constraints on the fastness 
parameter. For each source, we assumed M* = 1.4 M and = 1. We calculated j as a 
function of uu* from equation (33) and uf FM from equation (42) using the observed periods 
P*, period derivatives P*, X-ray luminosities Lx, and QPO frequencies z^qpo of these X-ray 
pulsars (Takeshima et al. 1994; Wang 1996, and references therein). We determined the 
mass accretion rates through M = LxR*/GM*. The magnetic dipole moments /i*3 have 
been recently deduced from the cyclotron features (Coburn et al. 2002). Sources like 4U 
1626 — 67, which are known to be close to their spin equilibrium, indicate lower values for 
the critical fastness parameter. The reason for this could be the deviation of angular velocity 
from the Keplerian profile near the inner disk regions (Li & Wang 1996). Figure 3a shows 
the predicted rotation curve of 4U 1626 — 67 when the pulsar was spinning up just before the 
torque reversal in 1990. After the torque reversal, we expect the angular motion of the disk 
matter around the same pulsar as in Figure 36. The estimated angular velocity profiles of 
4U 0115 + 63 and Cen X-3 are plotted in Figures 4 and 5, respectively. We could obtain the 
relevant rotation profiles shown in Figures 3, 4, and 5 only for specific values of oj* provided 
we have uj s = uf FM (see eq. [40]). With these constraints, the disk model depends on three 
parameters, r in , f3 and j, and there is a unique solution, such that if one of these three 
parameters is known, the other two can be obtained as the parameters that yield the unique 
solution for Q(r). For example, if the value of r in is assumed, e.g., r in = r\, then (3 and j, 
therefore the torque on the star can be obtained; numerically by scanning the /3,j parameter 
space until the unique solution is found. In the applications, we used observed values of Cl 
to determine j and scanned in (5 only, to use a shortcut to the model solution. We now 
summarize our results for each of these individual sources. 

4U 1626 — 67 is a low-mass X-ray binary pulsar with a 7.66 s spin period. The pulsar 
was in a steady spin-up state with P* = —4.86 x 10~ n s s _1 during 1977 — 1990. After the 
1990 torque reversal, the source underwent steady spin-down with P* = 4.32 x 10~ n s s _1 . 
Although there is no indication for an abrupt change in the X-ray luminosity at the torque 
reversal, a gradual decrease over a decade in the mass accretion rate can be seen from the 
archival flux history (Chakrabarty et al. 1997). The QPO frequency is ~ 0.04 Hz during 
both the spin-up and spin-down episodes. The estimated mass accretion rates are M 17 m 1 
for spin-up and Myj ~ 0.5 for spin-down. The beat frequency model gives wf FM = 0.76. 
Figure 3a shows the predicted rotation curve, when the pulsar spins up, for oo s = 0.76 with 
uo* = 0.65 and j3 = 6.48. The angular velocity profile during the spin-down episode was 
estimated for the same fastness parameter as in Figure 3b with u;* = 0.2 and (3 = 163. We 
find r in /r co = 0.75 and r in /r = 0.91 for spin-up, and r in /r co = 0.34 and r in /r = 0.41 for 
spin-down with r co = 6.5 x 10 8 cm. At the torque reversal, we expect that ^^s^ changes 
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from 4.7 to 3.8 according to our values of (3 and for /x^o = 2.2. Note that the angular 
velocity transition zone is quite broad for the spin-down state, i.e., 5r/r = 1 — r in /r = 0.59 
in comparison with 5r/r = 0.09 during the spin-up episode. 

4U 0115 + 63 is a 3.61 s pulsar. It has P* ^ -3.66 x 10~ 12 s s" 1 and M 17 « 0.43. The 
observed QPO frequency is 0.062 Hz. We inferred a fastness parameter of 0.817 from the 
beat frequency model. We find r = 3.4 x 10 8 cm and r co = 3.9 x 10 8 cm. The same value 
for uj s can be determined from the angular velocity profile if = 0.78 (see Fig. 4). We 
get r in /r = 0.97. The width of the transition zone can be estimated as 5r/r = 0.03. An 
Alfven radius of 4.3 x 10 8 cm for yU^o = 1 gives r in /rA = 0.77. The solution in Figure 4 was 
obtained for (3 = 4.25. This further implies that 'J^s^ = 1.7. 

Cen X-3 is a 4.83 s pulsar. It has P* = -4.29 x 10 -11 s s" 1 and M 17 « 4.3. The 
observed QPO frequency is 0.035 Hz. The fastness parameter inferred from equation (42) is 
0.855. We estimate r /r co = 0.9. A boundary layer width of 5r/r = 0.01 can be deduced 
from r in /r = 0.99. Our computation gives lj s = 0.855 only for uo* = 0.84 and (3 = 3.92 (see 
Fig. 5). This yields r in ^ 4.27 x 10 8 cm and ^ s 2 cS ^ 14 for ^ 30 = 2 and r co ^ 4.8 x 10 8 cm. 
Our model predicts r in /r A = 1.28 for this source. 

5. MODEL APPLICATIONS IN COMPARISON WITH RESULTS OF 

NUMERICAL SIMULATIONS 

The recent numerical study made by RUKL02 is the only work where the slow viscous 
accretion and the disk-magnetized star interaction were investigated in full two-dimensional 
time-dependent simulations with quasi-stationary conditions. The results of RUKL02 are 
relevant to all stationary models that describe the magnetic interaction of an accreting star 
with its surrounding disk. As the simulations by RUKL02 are two-dimensional whereas 
our work is limited to one-dimension, we expect significant differences between some of 
our results and their predictions. The inner boundary chosen by RUKL02 to investigate the 
dynamics of disk accretion is the rotating star on which free boundary conditions are applied 
to several hydrodynamic variables. In our work, the computational region is limited by the 
inner radius of the disk where the disk plasma corotates with the stellar magnetosphere. We 
keep the Keplerian law fixed as the asymptotic outer boundary condition on the rotation rate 
of the disk matter. The outer boundary conditions taken by RUKL02 on all hydrodynamic 
variables are fixed for the maximum computational region and free for a smaller simulation 
region (inner disk). 

Although our boundary conditions differ significantly from the ones used by RUKL02, 



-18- 



it is possible to make appropriate conversions between our reference values and those of 
RUKL02 for the inner disk radius and the physical quantities such as the stellar rotation 
rate and the angular momentum flux. Our reference value for the distance is always the true 
inner radius of the disk, r in , which is equivalent to one in dimensionless units. The innermost 
disk radius within which all the disk matter goes up into the funnel flow changes in time 
relative to the reference value for the distance, Rq, in numerical simulations of RUKL02. 
This innermost radius corresponds to r in in our model since the matter nearly corotates with 
the star inside r in < R . 

The dimensionless stellar rotation rate, that is the fastness parameter in our units is 
given by 



where is the fastness parameter in units of RUKL02. 

The dimensionless total angular momentum flux to the star in our units can be written 



reference mass accretion rate in units of RUKL02. 

In the following, we will estimate the approximate location of the braking radius r^ T , 
using our model parameters f3, j, and cj*. The braking radius was first defined by RUKL02 
as the radius inside which the disk is significantly disturbed by the stellar field. In almost 
all simulations by RUKL02 (especially the relevant ones with "Type I" initial conditions), 
the density in the disk drops sharply at r < r\„ and the matter is magnetically braked 
for r in < r < r br . The density in the disk rises again near the innermost radius. The 
radial variation of density at each time step of these numerical simulations is qualitatively 
reminiscent of the behavior of our dynamical viscosity, f(x) = 37rz/£/M (see Fig. 1). To 
make a quantitative comparison of our model predictions with the results of RUKL02, we use 
the numerical data of simulations for accretion to a slowly rotating star with = 0.19. We 
select the numerical values given by RUKL02 for the mass accretion rate M/M , the angular 
momentum flux tp, and the inner disk radius r m /Ro at times t = 30Po and t = AOPq, where 
Po is the rotational period at i?o- Using equations (43) and (44), we calculate our model 
parameters cu* and j. Equations (30) and (34) can be employed together with the condition 
f(xi) <C 1 (see §§ 3 and 4) to solve for (3. Once all model parameters are determined, the 
braking radius r br ~ (r in /R )xi where the density becomes minimum can be readily found 
from equation (34). Our numerical calculations yield r br ~ 2.61 at t — 30Pq and r br ~ 2.65 




(43) 



as 




(44) 



where J = ipM Rlfl K (R ) is the total angular momentum flux to the star and M is the 
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at t — 40P - These estimates are in agreement with the results of RUKL02. 

6. DISCUSSION 

We investigated the rotational dynamics of an accretion disk threaded by the dipolar 
magnetic field of a neutron star with constant screening factor. Within model assumptions, 
we found that a range of narrow or wide angular velocity transition zones with non-Keplerian 
rotation are consistent with observations. Our approach is analogous to the work of Glatzel 
(1992) for the non-magnetic case. We derived a vertically integrated dynamical viscosity form 
for a magnetically threaded disk from the conservation of mass and angular momentum using 
the Keplerian rotation as a particular solution for the angular motion of the disk plasma. 
Assuming that the whole disk can be treated by the same dynamical viscosity, we were able 
to obtain rotation curves without the restricting assumption of a thin boundary layer. We 
also used a more realistic, spatially varying azimuthal pitch, / B z , throughout the angular 
velocity transition zone instead of keeping it constant as in the model of Ghosh & Lamb. The 
effect of the magnetically enhanced viscous stress is included in our calculations. The rotation 
of the disk matter is almost Keplerian outside the corotation radius. The deviation from 
the Keplerian law is noticeable where the effective viscous stress is negligible as compared 
with the magnetic stress. This can be realized at a specific disk radius, r\ = £ir in , where 
the vertically integrated dynamical viscosity function, /, goes through a local minimum and 
actually nearly vanishes. The angular velocity at r\ is still nearly Keplerian. We find that 
the angular motion of the inner disk plasma from r\ to r in is controlled by both the viscous 
and magnetic stresses. The true inner disk radius, r in , is therefore determined in our model 
by the close balance between the magnetic stress and the viscous and material stresses. The 
width of the transition zone, 5 < x± — 1, depends basically on the model parameters f3, j, 
and uj* (see §§ 3 and 4). The physical meaning of this is that for different values of M, 
fi* and for given the torque on the star from a disk threaded by a screened dipole field 
can be achieved with the adjustment of the Q(r) profile through a transition zone. The 
rotation is nearly Keplerian throughout the whole disk and the angular velocity transition 
region is narrow in the case of dynamical viscosity decreasing towards the innermost edge of 
the disk for the current values of f3, j, and cj*. The boundary layer approximation of GL79 
remains valid only for f(x ~ 1) w (see Fig. 5). In this sense, the viscously unstable inner 
disk solutions of BC98 and CH98 could be an artefact of taking a fully Keplerian disk. In 
general the inner disk is not Keplerian; the angular momentum transfer with a dynamical 
viscosity increasing at small radii cannot satisfy the boundary condition imposed by the 
rotating magnetosphere if the inner disk is assumed to be Keplerian (see Fig. 36). 
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Using observed parameters of several sources in equilibrium, spin-up or spin-down states, 
we found that the non-Keplerian transition region for a Keplerian flow to match to the 
rotation rate of the stellar magnetosphere can be broad or narrow depending on source 
state. The beat frequency interpretation of QPOs provides a constraint on the fastness 
parameter for a given source, and thus on the model parameters (/3, j, a;*) that determine 
the appropriate rotation curve and the width of the transition zone. Among the examples 
we considered, the transition zones are narrowest, 5r/r ~ 0.01 — 0.1 for sources in spin-up 
(4U 1626 - 67, 4U 0115 + 63, Cen X-3); while 5r/r ~ 0.02 - 0.24 for sources in rotational 
equilibrium; and broadest, 5r/r ~ 0.6 for 4U 1626 — 67, when the source is in spin-down 
state. QPOs from disk-fed torque reversing X-ray binary pulsars may elucidate the function 
of non-Keplerian accretion flows and/or the width of the transition zones in the spin behavior 
of these sources. 

Results of the numerical simulations by RUKL02 were compared to our model. These 
two-dimensional simulations describe an analogous situation where the disk is threaded by 
the magnetic field of the star; both viscosity and diffusivity are of the same order of mag- 
nitude; and the disk matter slowly accretes across the field lines. Although the direct com- 
parison of the angular velocity distributions is not possible because of different conditions 
assumed at the boundaries of computational domains; there are remarkable similarities be- 
tween our model predictions and the results of RUKL02 (see § 5). 
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Fig. 1. — Radial variation of normalized vertically integrated dynamical viscosity f(x) = 
SnuT^/M from the outer disk regions x = 50 to the innermost radius x — 1, where x = r/r in , 
for different values of parameters (/3,j, a;*). The solid curve corresponds to f(x) with (3 = 16, 
cj* = 0.6, and j = —0.8 and the dashed curve to /(#) with j3 = 13, cj* = 0.3, and j = 0.9. 
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Fig. 2. — Radial variation of angular velocity u(x) = Q/Q^(r[ n ) and the corresponding 
normalized vertically integrated dynamical viscosity f(x) = 37rz/£/M for a disk around a 
magnetized star in spin equilibrium. The rotation rate profiles uu(x) shown by solid curves 
in panels (a) and (b) obtain for (3 = 14.5, oj* = 0.5, j = and (3 = 4.9, = 0.8, j = 0, 
respectively. In each panel the dashed curve corresponds to f(x) and the dotted curve to the 
Keplerian profile ujk(x) = x~ 3 ^ 2 , where x = r/r in . 
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Fig. 3. — Estimated rotation rate profiles u(x) = Q/Q^(r in ) and the corresponding normal- 
ized vertically integrated dynamical viscosities f(x) = 37rz/£/M for a magnetically threaded 
disk around an X-ray pulsar with observed parameters appropriate for 4U 1626 — 67 before 
and after the torque reversal. The panel (a) represents the angular velocity transition zone 
during the spin-up episode. The transition region is broad as seen from panel (b) when the 
source is in spin-down state. The rotation rate profiles cu(x) shown by solid curves in panels 
(a) and (b) obtain for (3 = 6.48, uo* = 0.65, j ~ 0.196 and (3 = 163, = 0.2, j ~ -0.513, 
respectively. In each panel the dashed curve corresponds to f(x) and the dotted curve to the 
Keplerian profile cok(x) = x~ 3 ^ 2 , where x = r/r- m . 
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Fig. 4. — Radial variation of angular velocity u(x) = fi/^K( r in) and the corresponding 
normalized vertically integrated dynamical viscosity f(x) = 37rz/£/M for a disk around an 
X-ray pulsar with observed parameters appropriate for 4U 0115 + 63. The rotation rate 
profile uo(x) shown by the solid curve obtains for (3 = 4.25, oj* = 0.78, j ~ 0.185. The dashed 
curve corresponds to f(x) and the dotted curve to the Keplerian profile ujk(x) = x~ 3 ^ 2 , where 
x = r/r in . 
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Fig. 5. — Radial variation of angular velocity u>(x) = f2/^K( r in) and the corresponding 
normalized vertically integrated dynamical viscosity f(x) = "invY^jM for a disk around an 
X-ray pulsar with observed parameters appropriate for Cen X-3. The rotation rate profile 
u>(x) shown by the solid curve obtains for (3 = 3.92, 00* = 0.84, j ~ 0.106. The dashed curve 
corresponds to f(x) and the dotted curve to the Keplerian profile u>k(x) = x~ 3 / 2 , where 
x = r/r in . 



